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ABSTRACT 

The role of HD cooling in the formation of primordial objects is examined by means 
of a great number of 1-D models of the collapse of halos, exploring a wide range of 
masses and virialization rcdshifts. We find that HD has very little effect upon the 
critical mass separating the objects which are likely to form stars from those which 
are not. We also find that, once the proto-stellar collapse has started, HD effects are 
quite negligible. 

Instead, HD effects can be important during the intermediate stage of gas fragmen- 
tation: objects below a certain mass scale (~3x 10 5 M Q at z v ; r = 20 in our "fiducial" 
case) can be cooled by HD down to T ~ 50 — 100 K, whereas H2 cooling never takes 
the gas below T ~ 200 K. The lower temperature implies a reduction of a factor ~ 10 
in the Jeans mass of the fragmenting gas, and stars forming in such low-mass halos 
are probably less massive than their counterparts in larger halos. We estimate the 
importance of this mode of star formation through a variation of the Press-Schechter 
formalism, and find that it never exceeds the contribution of halos which arc cooled 
by H2 only. Halos where HD is important account at best for a fraction ~0.25 of the 
total primordial star formation. However, HD cooling might provide a channel through 
which long-lived low mass stars could be formed in primordial conditions. 

Key words: stars: formation - molecular processes - cosmology: theory. 



1 INTRODUCTION 

Molecules play a decisive role in the formation of the first 
generations of luminous objects, because they provide the 
main cooling mechanism for primordial metal-free gas in 
small halos virializing at high redshift {e.g. Barkana & Loeb 
2001; Bromm & Larson 2004; Glover 2004; Ciardi & Ferrara 
2005; Ripamonti & Abel 2005). 

The H2 molecule has long been recognized (Saslaw & 
Zipoy 1967, Peebles & Dicke 1968, Matsuda, Sato & Takeda 
1969) as the most important cooling agent in such condi- 
tions. Its cooling and chemical properties have been carefully 
studied {e.g. Hollenbach & McKee 1979; Palla, Salpeter & 
Starrier 1983 [PSS83]; Stancil, Lepp & Dalgarno 1996,1998; 
Abel et al. 1997 [A97]; Galli & Palla 1998 [GP98], Glover, 
Savin & Jappsen 2006; Hirata & Padmanabhan 2006) and 
are a key component of all recent theoretical models and 
numerical investigations of the formation of the first ob- 
jects {e.g. Tegmark et al. 1997 [T97]; Omukai & Nishi 
1998; Bromm, Coppi & Larson 1999 [BCL99], 2002 [BCL02]; 
Abel, Bryan & Norman 2000, 2002 [ABN02] ; Ripamonti et 
al. 2002 [R02], Yoshida et al. 2006 [Y06]). 

Such studies suggest that primordial stars are more 
massive than their present-day counterparts, and that this 
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difference arises from the different cooling properties of 
metal-free and metal-enriched gas (Bromm et al. 2001; 
Omukai & Palla 2001, 2003; Schneider et al. 2002). A direct 
link seems to exist between the properties of H2 cooling and 
the large fragment mass and the high proto-stellar accretion 
rate which are found to occurr in simulations of primordial 
objects. 

However, most of these studies consider H2 as the only 
molecular coolant of primordial gas. This assumption might 
be quite critical, because of the close link between molecular 
and stellar properties. In fact, a number of chemical models 
of the primordial medium (Puy et al. 1993; GP98; Stancil, 
Lepp & Dalgarno 1998) have actually shown that, if the 
primordial gas cools below ~ 200 K, then HD molecules 
are likely to become the main cooling agents, despite their 
low number abundance. The cooling properties of HD and 
H2 are quite different from each other, so this could have 
important consequences: for instance, in an HD-cooled gas, 
the Jeans mass (and the typical mass of primordial objects) 
would likely be reduced. 

Such a possibly important role prompted further inves- 
tigations of the HD cooling and chemistry (Flower 2000; 
Galli & Palla 2002 [GP02]; Lipovka, Nunez-Lopez & Avila- 
Reese 2005) and especially of its effects on the properties 
of the first generations of stars (BCL99, BCL02; Uehara & 
Inutsuka 2000; Flower & Pineau des Forets 2001; Nakamura 
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& Umcmura 2002; Johnson & Bromm 2006; Lipovka et 
al. 2005; Nagakura & Omukai 2005; Shchekinov & Vasiliev 
2006; Y06). 

All the studies which find that HD effects might be 
important were considering situations where the primordial 
gas was already "perturbed" by some kind of phenomenon, 
such as the shock due to the merger of two halos, or the 
radiation from the first stars. The "unperturbed" case was 
considered by the detailed simulations by BCL99, BCL02, 
and Y06, who included a treatment of HD chemistry and 
cooling in some of their simulations, and found that HD 
had minor effects at best. However, such results concern the 
evolution of single "typical" cases and do not rule out the 
existence of objects where HD cooling is actually important. 

In this paper, we address the question of whether HD 
can be important during the birth of the first generation of 
luminous objects (that is, in unperturbed gas), with regard 
to the properties of the halos where they form, to the prop- 
erties of fragmentation inside such halos, and to the charac- 
teristics of the proto-stellar collapse. This is done by means 
of spherically symmetric ID calculations of the evolution of 
gas properties during the collapse of halos or proto-stars, 
through which we carry out an exploration of a large por- 
tion of the z v ir — Mhaio parameter space. Such exploration, 
combined with results from analytical models, allows us to 
give a quantitative estimate of the importance of HD-cooled 
objects. 

This paper is organized as follows: in section 2 we de- 
scribe the details of the ID calculations; in section 3 we 
examine the results obtained from such calculations, and in 
section 4 we discuss their cosmological significance. Finally, 
in section 5 we summarize our conclusions. 

Throughout all of this paper, we assume a flat ACDM 
cosmological model with parameters taken from the three 
years WMAP results reported by Spergel et al. 2006: = 
0.76, fi m = 0.24, tt b = 0.042, Q DM = n m -fi b = 0.198, h = 
0.73. Furthermore, p ~ 1.88xl0~ 29 fr 2 gcm~ 3 is the critical 
density of the universe, 



2 METHOD 

In order to assess the effects of HD in a wide range of cos- 
mological environments, a method similar to that of T97 is 
adopted: rather than simulating a single "typical" object in 
great detail (as was done by BCL99, BCL02 and Y06), the 
"coarse" behaviour of a plethora of halos is probed, explor- 
ing the Zvir — M ha i parameter space. 



2.1 The ID Code 

The main tool employed in this paper is a ID Lagrangian 
hydrodynamical spherically symmetric code, which follows 
the evolution of primordial gas inside and around DM halos, 
from the recombination epoch until well after the DM has 
virialized. Such code is described in R02 (see also Thoul & 
Weinberg 1995; and Omukai & Nishi 1998). Here we just 
outline the changes that were made in order to adapt it to 
our present purposes. 



Table 1. List of considered reactions 



# 


Reaction 








Reference 


1 


H + + e 


> 


H + 


7 


A AT /(111 

A97/2 

/"I T~> no / T T -1 n, 

GP98/H1 


2 


H + 7 




H+ 


+ e~ 


/~1 t~i n /Tin/) 

GF98/H2 


3 


H + e~ 


— * 


H- 


+ 7 


A97/7 


4 


H~ + 7 


— * 


H + 


e~ 


GP98/H4 


5 


H + H~ 




H 2 - 


f e" 


A97/8 


6 


H~ + H+ 


— * 


H + 


H 


A97/16 


7 


H + H+ 




H 2 + " 


r 7 


GP98/H8 


8 


H++ 7 


— * 


H + 


H+ 


GP98/H9 


9 


H++ H 


-* 


H 2 - 


f H+ 


A97/10 


10 


H 2 + H+ 


-* 


H 2 + - 


t- H 


GP98/H15 


11 


D+ + e~ 


-* 


D + 7 


A97/2 a 












GP98/Dl a 


12 


D + 7 


-* 


D+ 


+ e" 


GP98/D2 i ' 


13 


H+ + D 


-» 


H + 


D+ 


GP02/5 


14 


H + D+ 


-* 


H+ + D 


GP02/6 C 


15 


D+ + H 2 


— ► 


HD + H+ 


GP02/2 


16 


HD + H+ 




D+ + H 2 


GP02/4 


17 


H + e- 




H+ 


+ e~ + e~ 


A97/1 


18 


H + H 




H+ 


+ H + e~ 


PSS83/9 


19 


He + e~ 




He+ + e~ + e~ 


A97/3 


20 


Hc+ + e~ 




He - 


r 7 


A97/4 


21 


He+ + e~ 




He++ + e~ + e~ 


A97/5 


22 


He+++ e- 




Hc+ + 7 


A97/6 


23 


D + H 2 




HD + H 


GP02/l d 


24 


HD + H 




D + 


H2 


GP02/3 d 


25 


H + H + H 




H 2 


f H 


PSS83/4 


26 


H 2 + H 




H + 


H + H 


PSS83/5 


27 


H + H + H 2 




H 2 - 


f H 2 


PSS83/6 


28 


H 2 + H 2 




H + 


H + H 2 


PSS83/7 



The reference codes give the paper from which the reaction co- 
efficient was taken, and the number identifying each reaction in 
the corresponding paper. 

a We adopted the GP98 (case B) recombination rate before the 
turn-around redshift of each halo, when we switched to the A97 
(case A) rate. 

b The coefficients given by GP98 for the ionization of H (and 
D) were transformed into reaction rates we could promptly use 
through eq. (1) of Sasaki & Takahara 1993. 

c As noted by Johnson & Bromm 2006, the reaction coefficient 
given by GP02 contains a typo; therefore, we actually use the 
original rate of Savin (2002). 

d When extrapolated to low temperature ( 100 K), these re- 
action rates become extremely large; since this is an artifact of 
the adopted fitting forms, and there is no experimental results 
at T < 170 K, we assume that at T < 100 K the rate coeffi- 
cients of these reactions remain constant at the value we obtain 
for T = 100 K. 

2.1.1 Chemistry 

The chemical network of the original R02 code was revised 
and expanded. Now it follows the evolution of 12 species 
(the 9 "original" species, i.e. H, H + , H~, H 2 , Hj, He, He+, 
He ++ , and e~; plus D, D + , and HD). The considered re- 
actions are listed in Table 1. The network includes all the 
reactions for hydrogen and deuterium species which are in 
the minimal model of GP98, the collisional ionizations of hy- 
drogen and helium, and two other deuterium reactions de- 
scribed by GP02 (these reactions were not part of the GP98 
minimal model). Finally, there are the 3-body reactions of 
H 2 formation from PSS83. We stress that reactions 3 and 5 
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(13 and 15) are the main formation channel for H2 (HD), 
unless the density is so high (p > KT 16 gem" 3 ) that 3-body 
reactions (25 and 27) dominate. 

As in the original code, the non-equilibrium chemistry 
of the considered species is solved at each time step through 
an implicit difference scheme. 

It should be pointed out that, even if the rates given 
in Table 1 are fairly updated, their choice remains quite 
slippery: the uncertainties range from 10-20 per cent up to 
one order of magnitude in the worst cases (reactions 3, 4, 
5, and 6, including the main formation channel of H2 ). 
Glover et al. (2006) and Hirata & Padmanabhan (2006) 
find that the large uncertainties in these ill known reaction 
rates can have important effects on cosmological predictions; 
but also that the formation of the first protogalaxies from 
cold, neutral gas is relatively insensitive to the choice of 
the rates for these reactions. If so, the most problematic 
rate in our chemistry calculations is likely to be the one 
for reaction 15 (which is part of the main HD formation 
route): the value given by GP02 agrees within 20-30 per 
cent with that from Wang & Standi (2002), but there exist 
measurements differing by almost a factor of 2 from these 
theoretical estimates (see fig. 3 of GP02, and fig. 14 of Wang 
& Stancil 2002). We will try to estimate the effects of this 
uncertainty by running a dedicated set of simulations (see 
Section 2.2.2). 

It should also be noted that the rates of all the reactions 
which are caused by the interaction with a photon (2, 4, 8, 
and 12) only keep into account the effects of the CMB, while 
neglecting any extra radiation field which might be present 
(see Hirata & Padmanabhan 2006 for a possible effect of 
extra radiation). 



2.1.2 Cooling 

The original code included the treatment of cooling from 
H2 roto- vibrational lines (including their radiative transfer) 
and from the CIE (Collision-Induced Emission) continuum. 
This was subject to several changes: 

(i) The treatment of H2 cooling was simplified in order to 
make the code faster: for this reason, we used the H2 cooling 
rate from GP98, and treated the effects of line radiative 
transfer by using the methods devised in Ripamonti & Abel 
(2004) (in particular, of their equation 22); such methods 
are quite approximate, because the effects of line transfer are 
estimated only from the local density, whereas they depend 
also on temperature and velocity gradients. But they are 
adequate for the purposes of this paper. 

(ii) The treatment of CIE cooling was switched off, as 
we never approach the density regime where it is important 
(n > 10 14 cm" 3 ). 

(iii) We included the cooling due to HD molecules, adopt- 
ing the results of Lipovka et al. (2005) , which consider also 
the effect of vibrational lines, and are valid in a wide range 
of temperatures (40 K < T < 2 x 10 4 K) and densities 
(1 cm -3 < «h < 10 8 cm -3 ; this range can be easily extrap- 
olated to both higher and lower densities) 1 . 



1 We actually used the polynomial fit provided by Lipovka et 
al. (2005); the fit is quoted to be close to their original results for 



(iv) As we are concerned with a cosmological scenario, we 
included the cooling (or heating) from Compton scattering 
of CMB photons, and from Lyman a line cooling of atomic 
hydrogen. The adopted cooling rates per unit volume (from 
Rybicki & Lightman 1979 and Dalgarno & McCray 1972, 
respectively) are 

A C (T,T 7 ) ~ 1.0 x 10~ 37 n e T 4 (T 7 — T) ergs -1 cm -3 (1) 

A H (T) ~ 7.5 x HT 19 e- 118348/T n e n H ergs" 1 cnT 3 (2) 

where T, T 7 (~ 2.725(1 + z)), n e and 71h are the temper- 
atures of the gas and of the cosmic background radiation 
(both in K), and the number densities of free electrons and 
H atoms (both in cm -3 ), respectively. 

(v) When the temperature of the gas is of the same order 
as that of the cosmic background radiation, the cooling rates 
from H2 and HD can be substantially modified (see e.g. fig. 
9 of GP02), as atoms and molecules become heating agents 
when T < Tj. We account for the effects of the CMB upon 
the H2 and HD cooling rates by estimating the net cooling 
rates through an adaptation of eq. (6) of Puy et al. (1993) 2 



A not ,£r 2 (T, T 7 ) 

A ne t,HD,j(r, Tj) 



,A H2 (T)[1 
: A H d(T)[1- 



r H 2 (T T^") j 
^hd(tt-7T7)i 



(3) 
(4) 



where Th 2 — 512 K (Thd — 128 K) is the excitation tem- 
perature of H2 (HD). The CMB also affects the Lyman a 
cooling, but this is completely negligible, as A H (T) is ex- 
tremely small for T ~ T 7 . 



2.1.3 Dark Matter 

The gravitational effects of Dark Matter (DM) were intro- 
duced in the code. The DM halo was assumed to be always 
spherically symmetric and concentric with the simulated re- 
gion, whose central part represents the investigated halo. A 
DM mass of Mdm = Mhaio^DM/^m is assumed to be within 
a certain truncation radius R tr , inside which the DM density 
profile is a truncated isothermal sphere with a flat core of 
radius R ZO rc', outside the truncation radius, the DM density 
is assumed equal to the cosmological average. So, 



PdmM = 



Pcore 

Pc OIC (r/R c 



{ PoOdm(1 + z) 3 
where the core density p CO rc i 



if r < -Rcorc; 
if R 

core 

if r > R tT , 



= Mdm 



47T 



-Re, 



\ -fi-corc / 



(5) 



(6) 



as can be obtained by equating the DM mass within R tr to 
M DM . 

The truncation and core radii were assumed to evolve 
in time, mimicking the results for the evolution of a simple 
top-hat fluctuation (see e.g. Padmanabhan 1993), and are 
described by the following equations 



T > 100 K; however, it is reasonably accurate also for 40 K < 
T < 100 K. 

2 In the Puy et al. (1993) paper, the sign of the argument of the 
exponential is wrong because of a typografical error. 
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Rtr(z) = 



-Rcore(-2 ; ) — 



( 3 M DM \ 

\4t Pth( z ) ) 



1/3 



-Rvir 
Rvir 



tjz) 



Rtr(z) 
.TLvii 



t(z vir )-t(z t a) 



(g-|)*M 



if 2 > 2 ta 

if Z t a > 2 > «v 
if 2 < Z v i r 



(7) 



t( Zvir )-t( Zta ) 



if 2 > 2 ta 

if 2 ta > 2 > 2 vir (8) 
if 2 < 2 v i r 



where 2ta — 1.5(1 + 2 v ir) — 1 is the turn-around redshift, t(z) 
is the time corresponding to redshift z, and the DM density 
inside the halo before z ta is given by equations (23) 3 and 
(24) of T97, which are a good approximation of the top hat 
results: 



pTn{z) 

A 



Po^dm(1 + . 



1 + e 1-0.75^ - 1 



A(Z) = (1 + Zvir)/(1+Z). 



(9) 
(10) 



Finally, £ (for which we will consider values in the range 
0.01 — 0.2) is the ratio between the final size of the flat 
density core and the virial radius, which we define as 



R 



= — Rt 



(Zta 



47T pTH(2ta 



1/3 



(11) 



Such a definition is unusual, but the difference with the most 
common definitions of the virial radius (e.g. Padamanbhan 
1993) is less than 0.4 per cent. 

The above equations divide the evolution of the DM 
profile in three stages: before the halo turn-around, it is 
reasonable to assume that the density profile within the halo 
is flat (therefore, Ticoro = Rtr). After the turn-around the 
density profile is evolved smoothly, reaching an equilibrium 
configuration at the virialization redshift; after that, the DM 
profile becomes completely static. 

This is a crude model of the dark matter halo, but it is 
adequate for our purposes. We emphasize that the choice of 
a model with a flat core, rather than a cusp as in the NFW 
profile (Navarro, Frenk & White 1997), helps to ensure that 
the behaviour we observe near the centre is due to the self- 
gravity and hydrodynamics of the simulated gas, rather than 
to the assumed DM profile. However, in the following we will 
also discuss the effects of a cuspy NFW profile for the static 
post-virialization phase. 



2.2 The models 

2.2.1 Initial conditions 

We start our computations at a very high redshift (2 = 
1000), when both the dark matter and baryonic overden- 
sities are still very small, and it is perfectly reasonable to 
assume an uniform distribution of the gas (the effects of 
Compton drag further improve the accuracy of this assump- 
tion - see eq. 4.171 of Padmanabhan 1993). Furthermore, we 
simulate a region whose comoving initial radius is 10 times 
larger than that of the halo, in order to achieve a better 
modeling of the hydrodynamical effects, as this guarantees 
that the total mass in the simulated region is well above both 

3 In the original reference (T97) this formula is affected by a 
typografical error, as the sign inside the exponential is wrong. 



Table 2. Assumed chemical composition at z=1000 



Spccies(X) 


"x/"H ° 


Comments 


H° 


0.9328 


< 1 because H is partly ionized 


H+ 


0.0672 


Sasaki & Takahara 1993 6 


H" 


1Q -19 


GP98, fig. 4 


H 2 


10- 13 


GP98, fig. 4 


H 2 + 


io- ls 


GP98, fig. 4 


He 


0.0833 




He+ 


10 -25 


GP98, fig. 4 


He++ 







e~ 


0.0672 


from charge conservation 


D° 


2.332 x 10" B 


(Romano ct al. 2003) c 


D+ 


1.68 x 10~ 6 


(Romano ct al. 2003) c 


HD 


2.5 x 10- 18 


(Romano ct al. 2003) c 



n H + ) + "HD ■ 
b = 0.05, h = 0.5 (tab. 1); this 



a iih is the total abundance of H, including all species 
™H = «H0 + n H+ + n H- + 2 ( n H 
6 from their model with Qo = 1, fi 
value also agrees with the result of the RECFAST code (Seager, 
Sasselov & Scott 1999, 2000) for the flat ACDM cosmology we 
are assuming. 

c Romano et al. (2003) actually give the total abundance of 
D, rt D = 2.5 X 10~ 5 n H ; we split this into the D°, D+ and HD 
abundances by assuming that n D + /tl-q = n^/tin, and that 
«HD/nD = «H 2 / n H- 



the cosmological Jeans mass and the filtering mass (see e.g. 
Peebles 1993, and Gnedin 2000). 

We assumed that at z — 1000 the gas temperature 
is the same as that of the CMB (~ 2728 K), that its 
density is equal to the cosmological average for baryons 
(p gas (2 = 1000) ~ 4.22 x 10~ 22 gem" 2 ), and that the chem- 
ical abundances are those listed in Table 2. 



2.2.2 Simulation sets 

We ran several sets of simulations (listed in Table 3) , cover- 
ing halos with virialization redshifts in the range 10 < 2 v ir < 
100 and total (baryonic+dark matter) halo masses in the 
range 5 x 1O 3 M < M ha io < M H (2 v ir). M H is the minimum 
mass of the halos in which the cooling from atomic H is effec- 
tive, that is, the minimum mass of halos with T v i r > 10 4 K 



M H (2 vir ) ~ 1.05 x 10 9 M C . 



(JL.) 

V 1.237 



-3/2 



(1 + 2 v ir) 



-3/2 



(12) 



where fi is the mean molecular weight of the primordial gas. 

Each simulation set includes two runs for each 
(2 v ir, Mhaio) pair, differing only because HD cooling is ei- 
ther included or omitted. The results of the runs without HD 
cooling are used as a control sample. The models are based 
on the initial conditions described above, and consist of 150 
shells, whose mass increases from the centre to the outskirts. 
The spacing was chosen in order to achieve sufficient mass 
"resolution" at the centre: the mass of the central shell is 
always ~ 0.3M Q , which is a very small fraction (~ 10~ 3 ) of 
the mass of the fragments seen in the simulations of ABN02, 
BCL02, and Y06. 

We ran several sets of simulations, in order to check 
the effects of the most uncertain parameters describing our 
initial conditions, such as the DM density profile: 
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Table 3. Summary of simulation sets 



Set 


DM profile 






Abundances 


Remarks 


£ = 0.1 


Isothermal" 


« = 


0.1) 


from Table 2 


Fiducial 


cxtra-D 


Isothermal" 




0.1) 


High c D, e", H+ 


Doubled rate of reaction 15 


£ = 0.2 


Isothermal™ 




0.2) 


from Table 2 




£ = 0.05 


Isothermal™ 


« = 


0.05) 


from Table 2 




£ = 0.01 


Isothermal" 


(? = 


0.01) 


from Table 2 




NFW 


NFW 6 (c = 


10) 




from Table 2 





a Truncated isothermal, as described in section 2.1.3. 

b Before virialization, the DM profile is assumed to evolve as in the fiducial (§ = 0.1) set. 

c For D, e — , and H + , twice the value given in Table 2. The other abundances given there are the same, except for H° (which was 
reduced to 0.8656 in order to compensate the increase in H+). 



(i) We check the influence of the concentration of the 
previously described truncated isothermal DM density pro- 
file by running four sets (£ = 0.2, £ = 0.1, £ = 0.05, and 
£ = 0.01) of simulations differing for the value of the £ pa- 
rameter. We will refer to the £ = 0.1 set as to the "fiducial" 
set. 

(ii) We check whether a cuspy density profile can lead 
to significantly different results by running a set of simula- 
tions (NFW) in which the post-virialization density profile 
within _R v i r is described by a NFW provile with concentra- 
tion parameter c = 10; in these models, the evolution of the 
DM density profile before virialization is the same as for the 
"fiducial" model. 

(iii) We check whether HD effects depend strongly on the 
assumed abundances (at z = 1000) and on the uncertain rate 
of reaction 15 by running a set of simulations (extra-D) in 
which we double the initial values of both the electron (and 
H + ) fraction given in Table 2 (to 0.1344), and the deuterium 
total abundance (to nu/nn — 5 x 10 -5 ); we also assume that 
the rate of reaction 15 is twice as high as its standard value. 
All these changes favour HD formation, and this simulation 
set is mostly useful as an upper limit on the effects of HD 
(in the case that £ = 0.1). 

All the runs were stopped when a density p co n = 
1.67 x 10 -19 gem -3 = 10 5 mH cm -3 is reached at some shell 
(generally at the centre), or when the simulation reaches 
Zstop = 5, whichever comes first. Our choice of p co n ensures 
that the gas reaches densities where the first run-away col- 
lapse is taking place; furthermore, p co n 3> p v ir, which en- 
sures that the collapse is highly significant, and also that 
our results are quite independent of the exact value of p co ri. 
Instead, we chose z st0 p = 5 because below that redshift the 
halos which are cooled by molecules are both cosmologically 
unimportant, and extremely unlikely to survive in the un- 
perturbed (and neutral) conditions we are assuming. 




1000 500 300 200 100 50 30 20 

Redshift 



Figure 1. Chemical and thermal evolution between recombina- 
tion and virialization of a IO^Mq halo with z vir = 20, taken from 
the "fiducial" set of simulations. The left axis refers to the labelled 
lines, showing the abundance evolution of the various species: the 
abundances of H+ and e~ are always extremely close, and are 
represented by a single line; neutral and ionized He were omitted. 
The right axis refers to the two dot-dashed lines without labels, 
which show the evolution of the CMB and of the gas temperature 
(thin line and thick line, respectively). The results shown in this 
plot are almost independent of Mhaio and of the inclusion of HD 
cooling, up to at least the turn around redshift Zt H = 30.5. 



which will later virialize into halos, such as those shown in 
Fig. I 4 . 

Instead, HD cooling might affect the thermal evolution 
of the gas in a halo before and immediately after virializa- 
tion, changing the critical mass which separates efficiently 
and inefficiently cooling halos (see e.g. T97). 

Each of our simulated halos was classified as collapsing 
(or equivalently, efficiently cooling) or non- collapsing (inef- 



3 RESULTS 

3.1 Critical masses 

Since the Compton heating due to residual free electrons 
largely dominates the thermal behaviour of the "average" 
IGM, the inclusion of HD cooling scarcely affects the chem- 
ical and thermal evolution of the IGM, even inside regions 



4 Hirata & Padmanabhan 2006 have shown that the amount of 
H2 which formed in the IGM through the H^ and the H — chan- 
nels (reactions 3, 5, 9, and 10 of Table 1) at z ^ 70 might be 
substantially less than what is shown in Fig. 1. In the rest of this 
paper we neglect this uncertainty because in our simulations the 
bulk of H2 forms inside the halos after they virialize, and the 
"background" abundance has very little effects on our results. 
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Figure 2. Evolution of the critical mass (M cr i t ), and of the mass 
below which HD becomes important (AfjjD) as a function of the 
virialization redshift z v - lT , for our six sets of simulations. Solid lines 
show the evolution of M cr i t when HD cooling is included in the 
simulations, while dashed lines show the case where HD cooling 
is neglected; dotted lines show Mhd- Panels refer to the (fiducial) 
£ = 0.1 set (a), to the extra-D set (b), to the £ = 0.2 set (c), to 
the § = 0.05 set (d), to the § = 0.01 set (e), and to the NFW 
set (f). The results of the fiducial case are repeated as thin lines 
in all the panels, in order to facilitate the comparison between 
the various sets. M cr ; t is almost unaffected by the inclusion of 
HD cooling; however, the mass range between A/jjd an d M cr ; t 
becomes significant in the extra-D set, or if the DM halos are 
quite concentrated. 

ficiently cooling), depending on whether it reaches a maxi- 
mum density larger than p co n in a less than an Hubble time, 
i.e. at a redshift z co n ^ [0.63(1 + jz v ir)] — 1. 

Efficiently cooling halos have a high probability of form- 
ing luminous objects; this probability is much lower (but 
sometimes not completely negligible, as we show in the next 
section) for inefficiently cooling halos. 

Our definition of collapsing halos is comparable to the 
"collapse criterion" of T97 (see their §5.1), so we use a sim- 
ilar notation and define M cr it(z v ir) as the minimum mass of 
an efficiently cooling halo virializing at z v i r . 

In figure 2 we show the behaviour of Af cr it(zvir), com- 
paring the results of the runs including HD cooling with 
those of the "control" runs. All the six sets of simulations 
are shown. In all cases HD has very little effect upon M cr it, 
as the solid and dashed curves of each simulation set are 
always close: models where the DM is quite concentrated 
exhibit slighly larger differences; but they never amount to 
more than 20 — 30%, even in the case of a NFW profile. 

The reason of this result is that HD cooling can over- 
come H2 cooling only when the temperature is quite low 
( 200 K, see fig. 3). But the gas temperature at virializa- 
tion is T ~ T v i r ~ 1000 K, and HD can become important 
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Figure 3. Ratio of the HD cooling per molecule to the H2 cooling 
per molecule as a function of temperature at three different densi- 
ties: 1 cm -3 (i.e. low density limit; solid line), 10 4 cm" 3 (dashed 
line), and 10 s cm -3 (i.e. high density limit; dot-dashed line). 
They are compared to the njj 2 /iHD ratio (dotted line) which 
can be expected if reactions 15 and 16 of Table 1 are at equilib- 
rium (see e.g. of Shchekinov & Vasiliev 2006): if this assumption 
is correct, the intersections between this line and the other three 
lines mark the conditions where H2 and HD are equally important 
coolants. 

only after H2 has lowered the gas temperature by a substan- 
tial factor, i.e. only if H2 cooling is efficient in the first place. 
Such an explanation is confirmed by Figures 4 and 5: before 
virialization HD cooling might even exceed H2 cooling, but 
remains negligible when compared to the effects of Compton 
scattering on free electrons; at virialization H2 overcomes 
both the Compton (because the increase in density favours 
H2 formation, reducing the number of free electrons at the 
same time) and the HD contribution (because the increase 
in temperature reduces the ratio Ahd /Ah 2 , an d also the ra- 
tio fiHD/nH 2 )- So, HD cooling may become dominant only 
well after virialization, when H2 has already reduced the 
temperature by a factor 2. 

3.2 Fragmentation 

The simulations clearly show that even if HD has very little 
effect on which halos cool, collapse, and form stars, it can 
affect how this happens. Actually, there are objects where 
HD cooling has no effect, and objects where the thermal 
evolution of the gas is deeply affected by HD. 

In Figures 4 and 5 we show the evolution of the central 
properties in two halos chosen from the fiducial set. These 
halos have the same z v i r (20) but different mass (Mhaio = 
2 x 10 5 , 1O 6 M ). In each fi gure the results of runs with 
and without HD cooling are compared. In particular, we are 
interested in the evolution of the Jeans mass, which is taken 
to be (cfr. Peebles 1993) 

, T \ k ( 7rfc B T V /2 _i/2 

~ 5OM T 3/2 /i' 3/2 n- 1/2 . (13) 
In the case of the smaller halo (fig. 4) the final stages 
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Redshifl Redshift 

Figure 4. Evolution of the properties of the gas at the centre of 
an halo with M^aio = 2 X 1O 5 M0 virializing at z v i r = 20 (taken 
from the fiducial set), as a function of redshift. In all panels, thick 
lines refer to models where HD cooling was included, whereas 
thin lines refer to models neglecting HD cooling. The properties 
shown are: Jeans mass (panel a); gas density and temperature 
(panel b; solid and dashed lines, respectively); fraction of the total 
cooling/heating due to HD, H2 and Compton scattering (panel c; 
solid, dashed and dotted lines, respectively); fractional abundance 
of HD, H2 and electrons (panel d; solid, dashed and dotted lines, 
respectively; the HD abundance shown here is 10 3 times larger 
the real value). In panel (b) the value of the CMB temperature is 
also shown (dotted line), and in panel (d) we additionally show 
the evolution of the jihd/'I'H^ ratio (dot-dashed line). 

of evolution differ, because HD cooling becomes dominant, 
anticipating the collapse (in this case p co n is reached at 
z ~ 11.5 rather than z ~ 10.3), and, most importantly, 
lowering the gas temperature. In fact, the final temperature 
is ~ 70 K, i.e. 3-4 times lower than in the case where HD 
was neglected; as the Jeans mass is reduced by a factor ~ 10, 
this behaviour strongly hints towards a change in the frag- 
mentation properties of the gas inside an halo of this kind. 
On the other hand, in the case of the larger halo (fig. 5), the 
presence of HD makes no difference. 

Figures 4 and 5 are illustrative of a general trend: HD 
is unimportant in the largest halos, but it becomes relevant 
at lower masses. In fact, when the largest halos are consid- 
ered, the H2 abundance just after virialization is above the 
threshold (~ 5 x 10~ 4 ) identified by previous studies (e.g. 
T97) for H2 cooling to be efficient. As H2 sheds away most of 
the compressional heating, the collapse proceeds unimpeded; 
the gas remains at temperatures T above the 200 K "thresh- 
old" below which HD cooling becomes important (T usually 
reaches a minimum between 200 and 300 K, then slowly in- 
creases as the collapse proceeds). Instead, in smaller halos 
the H2 abundance is below the threshold for efficiently dis- 
persing the compressional heating, and the gas undergoes a 
stage of slow contraction and cooling; for example, in fig. 4 
this phase lasts from z ~ 20 to z ~ 13. The "extra" time 




26 25 24 23 22 21 20 19 26 25 24 23 22 21 20 19 

Redshifl Redshift 

Figure 5. Evolution of the properties of the gas at the centre of 
an halo with Mhaio — 10 6 Mq virializing at z v i r = 20, (taken from 
the fiducial set) as a function of redshift. See the caption of fig. 4 
for a description of the panels and an explanation of the symbols. 
The absence of the thin lines associated with the case where HD 
cooling was not included is only apparent: in halos as massive 
as the one shown here, the gas evolution is not influenced by HD 
cooling, and the thin and thick tracks are perfectly superimposed. 



spent in this phase allows a build up of HD (see especially 
panel (d) of fig. 4), until HD completely dominates the cool- 
ing and is finally able to "restart" the collapse. 

In figure 2 we compare the critical mass M cr i t (z) to the 
mass Mhd(z) below which HD cooling becomes important 
(Mhd(z) was taken as the maximum mass where the inclu- 
sion of HD cooling takes the final gas temperature below 
200 K, with a reduction of at least 25 per cent when com- 
pared to the value obtained when HD is neglected), for all 
the simulation sets we considered. The range between Mhd 
and M cr it (i.e. the collapsing halos where HD is likely to 
be important) is negligibly narrow for £ = 0.2, but it gets 
wider when an higher D abundance is considered, or when 
the DM halos are assumed to be more concentrated than in 
the fiducial case (see also Section 4). 

It is important to remark that our results are not in con- 
tradiction with the detailed simulations of BCL99, BCL02, 
and Y06, because all of them considered halos with masses 
in the range 0.5 — 1 x 10 6 Mq at redshift ~ 20, i.e. above 
Mhd- We must also remark that our results depend on the 
assumption that halos remain unperturbed, because mergers 
lead to dynamical heating, which could prevent the collapse 
also in halos with mass larger than M cr i t (see Yoshida et 
al. 2003; Reed et al. 2005). Since the probability of mergers 
is particularly high in halos where HD is important (because 
of the long "build-up" phase), this effect is likely to reduce 
the fraction of halos where HD cooling is important (see also 
Section 4.2). 
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3.3 Proto-stellar collapse 

Lipovka et al. (2005) suggested that HD cooling might be 
relevant also at relatively high densities and temperatures 
(T ;$ 3000 K, n ^ 10 6 ), a regime which is associated with 
proto-stellar collapse. Another reason to check this regime is 
that HD cooling might be in the optically thin regime even 
when the H2 cooling is substantially reduced by radiative 
transfer effects. 

In order to evaluate these effects, we re-ran some of 
our models with a better resolution (200 shells, mass of the 
central shell ~ 10 -3 Mq) stopping them only after a density 
n = 10 13 cm -3 was reached (after that, the fast increase of 
CIE cooling will overwhelm both H2 and HD line cooling, 
see e.g. Ripamonti & Abel 2004). We explored 6 virialization 
redshifts (z v ir = 10, 15, 20, 30, 50, 100), always choosing 
an halo mass of 10 6 Mq, above both M CI - lt {z) and Mud(z)): 
here we ignore halos where HD can affect fragmentation, 
because otherwise we cannot have a meaningful comparison 
with a model where HD cooling was neglected. 

Again, by comparing simulations with and without HD 
cooling, (see e.g. fig. 6) we find that HD causes no differ- 
ence also in this phase of the formation of a primordial 
star. The main reason is that the ratio jihd/ih2 decreases 

" -465K 

with temperature (at equilibrium, nHD/fiH 2 oc e t , see 
Shchekinov & Vasiliev 2006), and so does the ratio of 
HD cooling per molecule to H2 cooling per molecule. Fur- 
thermore, radiative transfer effects become important only 
when the H2 fraction is already of the order of 1, implying 
hhd/«h 2 ~ 2riD/fiH — 5 x 10 -5 (10 -4 if the abundances 
of the "extra-D" set are considered): this difference is sim- 
ply too large to be compensated by the radiative transfer 
effects on H2 cooling efficiency, especially at the relatively 
high temperatures (~ 1000 K) typically associated with the 
high-density phases of the collapse. 

We emphasize that the above calculations always as- 
sumed the optically thin limit for HD cooling. We have 
not checked whether this is correct, but our results will not 
change even if it is not: in fact, the treatment of radiative 
transfer would then yield an effective HD cooling rate which 
is lower than what we assumed. 



4 DISCUSSION 

In the previous section we showed that in the standard 
ACDM model the main effect that HD might exert upon 
the evolution of unperturbed primordial halos is to re- 
duce the typical stellar mass in objects below a certain 
mass MHD(zvir). Here, we investigate whether this "HD- 
dominated" regime of star formation is cosmologically rele- 
vant or not. This is done by comparing the mass of the stars 
which formed in halos affected by HD cooling, and in halos 
where only H2 (or H) cooling was important. 



4.1 Rates of formation and survival of objects 

We base our estimates on a variation of the Press-Schechter 
formalism (Press, Schechter 1974), namely the extended 
Press-Schechter (EPS) formulae of Sheth & Tormen (1999) 
[ST99] , which was found to be in reasonable agreement with 
the results of numerical simulations of the formation of very 
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Figure 6. Evolution of the properties of the gas at the centre of 
an halo with M^aio — 1O 6 M0 virializing at z v i r = 30 (taken from 
the fiducial set), as a function of central density. Panels (c) and 
(d) are the same as in figs. 4 and 5, except that the X-axis shows 
the gas central density, rather than the redshift. Panel (a) shows 
the time at which each density is reached, expressed as the time 
(tl3 — t) remaining before reaching a density 10 13 ran cm -3 at 
time ti3. Panel (b) shows the temperature evolution of the gas 
(solid lines) and of the CMB (dotted line). As in fig. 5, the absence 
of the thin lines (corresponding to the case where HD cooling was 
not considered) is only apparent: the results with and without HD 
are almost identical. Finally, the increase in the fraction of total 
cooling due to HD at high densities (panel c) is due to radiative 
transfer effects reducing the H2 cooling efficiency; such effects are 
not considered for HD, so the high density part of the solid curve 
in panel (c) should be considered an upper limit. 



early objects by Gao et al. (2005; see also Jenkins et 
al. 2001, Reed et al. 2003, Gao et al. 2004, Springel et 
al. 2005). 

Such formulae can be treated in the same way as Sasaki 
(1994) did for the original Press-Schechter, leading to similar 
results: the comoving number density of virialized halos in 
the mass range between m and m + dm at redshift z is 



N ST (m, z) 



So 




m 





it a 



a(m) 2 dm 



AS C 



x ^7-^ — e 2 

D(z) 



(14) 



where a = 0.707, A — 0.322 and p = 0.3 are the parameters 
given by ST99 (instead, a — 1, A — 0.5, p — correspond to 
a "standard" Press Schechter function), a(m) is the value of 
the root mean squared fluctuation in spheres that on average 
contain a mass m, S c ~ 1.686 is the overdensity threshold 
for the collapse, D(z) is the growth function of perturbations 
(see Peebles 1993), and 



(15) 



D(z) 2 a(m) 2 
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The time derivative of eq. (14) is 

N ST (m,z) = 



dDdz 1 
~dz~dt~D 

x NsT(m, z 



1 — av — 2p 



1 + (au)-p 



(16) 



and by following Sasaki 1994 (in particular in the assump- 
tion that the destruction rate of halos has no characteristic 
mass scale) we get that the formation rate is given by 



■Nform (m, z) 



dD dz 1 av + (av 



1-P 



2p 



dz dt D 



1 + (av)~P 



N ST (m,z) (17) 



and also that the probability that an object which exists at 
z survives until redshift z (z < z') without merging is 



p(z',z) = 



\D(z') 


l-2p 


n + z~ 


l-2p 


[D(z)_ 




ll + z'_ 





(18) 



where the last equality is strictly valid only in an Einstein-de 
Sitter universe with S7 m = 1. 



4.2 Mass fractions 

The results of our simulations can be interpolated in order to 
get an estimate of the function z co n(m, z v j r ), i.e. of the red- 
shift at which the collapse of the gas inside an halo of mass 
m that virialized at z v ir has proceeded far enough for star 
formation to occur; from the same data, we can also estimate 
a correlated quantity, i.e. the minimum mass M m i n (z v ir, z) 
that an unperturbed halo which virialized at z v ir must have 
in order to have formed stars before redshift z. These func- 
tions can be combined with the above equations (17) and 
(18) in order to find the mass fractions in object-forming 
halos of various kinds (HD-cooled, HVcooled, H-cooled) as 
a function of redshift. For example, the fraction of mass that 
at redshift z is in halos which have produced (or are produc- 
ing) stars through HD-dominated cooling is 



/hd(«) 



dz V u 



dt 



where Z\ = 100 is the maximum considered redshift (we 
assume that the contribution to the star-forming density of 
halos from redshift higher than z\ is negligible), and 



M H d,i = min(M H D(zvir), M min (z vir , z)), 

M H D,2 = M H D(Zvir) 



(20) 
(21) 



We note that the p(z v i T ,z) factor in eq. (19) excludes 
all the halos which experienced a major merger between 
Zvir and z; this is a very approximate way of accounting 
for the effects of dynamical heating (Yoshida et al. 2003; 
Reed et al. 2005) . Some of these excluded halos (namely, the 
ones undergoing a merger at a redshift between their z co ii 
and z) actually formed stars in the HD regime. This effect 
could be included in the calculation by changing the proba- 
bility factor into p(z v ir, max(z, z co ii(m, Zvir)))- However, the 
difference between the two probabilities is relatively small: 
in the most extreme case (i.e. z v ir = z co ii = 100; z = 10) 
p(z V ir, z co n)/p(z vir , z) ~ 2.41, whereas for much more typi- 
cal values of z v i r , z coU , and z the correction factor is ^ 1.5, 
which is too small to alter the qualitative conclusions that 
can be reached through eq. (19). 
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Figure 7. Evolution of the fraction of mass inside halos which 
were able to form stars within a given redshift. Solid lines: mass 
inside H2 -cooled halos; dashed lines: mass inside HD-cooled ha- 
los; dotted lines: mass inside H-cooled halos; dot-dashed lines: 
ratio between the masses inside HD-cooled halos and H2 -cooled 
halos. Panels refer to different sets of simulations, and are in the 
same order as in Fig. 2. 



The mass fractions /h 2 ( 2 ) an d /h(z) (due respectively 
to objects where H2 cooling and atomic H cooling is dom- 
inant) can be found by simply changing the limits of the 
mass integration inside eq. (19) to 



dm — A/f orm (m, z vir )p(z v ir, z) 
Po 



(19) 



Mh 5i1 = 

Mh 2 ,2 ~- 

and to 



min(M H (z vir ), max(Af H D(2vir), M min (z vir , z))X,22) 



M H (z v ir) 



M H ,i = max(M H (z v ir), M mill (z vir , z)), 

M H ,2 = OO. 



(23) 

(24) 
(25) 



We calculated /hd, /h 2 an d /h by taking the er(m) 
and D(t) functions given in Eisenstein & Hu 1998. In fig. 7 
we compare their redshift evolution. It is apparent that 
/hd < /h 2 a t a-U redshifts and in all the cases we are con- 
sidering. However, there is a substantial difference between 
the various simulation sets. In the least concentrated case 
(£ = 0.2) the amount of star formation through the HD 
cooling regime is likely negligible, never exceeding 4 per cent 
of the total and rapidly decreasing when going to high red- 
shifts. Instead, when the DM profiles are assumed to be 
moderately or highly concentrated (£ ;$ 0.05, or NFW pro- 
file), the difference between the HD and the H2 channels for 
star formation reduces to a factor ~ 3; it is also essentially 
constant between z = 10 and z = 30, even if the efficiency of 
the HD channel rapidly drops at z ^ 30 (a behaviour which 
could be expected from fig. 2); the fiducial set is somewhat 
in the middle, while the results from the extra-D set resem- 
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ble the ones we obtained for moderately concentrated DM 
profiles, except that the drop in /hd / /h 2 is much slower. 



5 SUMMARY AND CONCLUSIONS 

In this paper we examined the possible effects of HD cooling 
upon the most popular scenarios for the formation of first 
objects. 

We found that HD cooling has very little effect upon the 
evolution of a halo before and during virialization, so that 
the critical mass M cr i t separating star-forming and non-star- 
forming halos is scarcely affected. Similarly, we found that 
HD is even less important when the first phases of proto- 
stellar collapses are considered. Both these conclusions are 
quite independent of the shape of DM profiles, and of the 
exact abundance of Deuterium. 

The most interesting of our results is that HD can in- 
fluence the fragmentation process of primordial gas clouds 
into stars. In fact, our simulations show that in halos with 
masses just above M cr i t , HD cooling dominates the phases 
of gas collapse when fragmentation is likely to take place, 
significantly reducing the gas temperature, the Jeans mass, 
and probably also the typical mass of fragments. If so, HD 
cooling could open a channel for the formation of relatively 
low mass stars even in metal-free gas. 

We then employed a simple analytical model, based on 
the extended Press-Schechter formalism, in order to estimate 
the importance of this "HD mode" of star formation. A com- 
parison with the common "H2-only" mode of primordial star 
formation shows that the HD mode is always sub-dominant. 
However, its importance depends on the detailed DM den- 
sity profile, on D abundance, and on the rate of HD forma- 
tion. If DM profiles are assumed to be relatively loose, and if 
we take the "standard" values for D abundance and HD for- 
mation rate (as in our simulation sets £ = 0.2 and £ = 0.1), 
the range of halo masses where HD is important is quite nar- 
row, and the number of stars whose formation is triggered 
by HD cooling is probably negligible. Instead, in the case of 
quite concentrated profiles, and/or of higher D abundances 
and HD formation rate (such as in our £ = 0.05, £ = 0.01, 
NFW, and extra-D simulations sets), HD dominated halos 
should account for about one quarter of all primordial star 
formation, at least if we are not underestimating too much 
the effects of dynamical heating. 

We remark that there exist at least two mechanisms 
which could boost the importance of HD for primordial star 
formation. First of all, it is not unreasonable to expect that 
the star formation efficiency is higher for the halos which 
form stars through HD than for those where HD is unim- 
portant. This is because the lower typical stellar mass should 
imply a weaker feedback. 

A second possibility comes from considering that the 
ratio of HD to H2 (and with it the mass range where HD 
cooling is important) is quite enhanced when the abundance 
of free electrons is higher than in the standard case (see e.g. 
Nagakura & Omukai 2005, O'Shea et al. 2005, Johnson & 
Bromm 2006, Yoshida 2006, even if all of them consider a 
different scenario). This might happen because of DM an- 
nihilations or decays (Ripamonti, Mapelli & Ferrara 2006), 
or, less exotically, because of the influence of nearby ioniz- 



ing sources (e.g. accreting black holes, as in Zaroubi & Silk 
2005). 
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